
### MONTE CARLO SIMULATION FOR HIERARCHICAL MODEL
###
### /code folder individually has separate code for each simulation
##  "G=10_J=100_T=40.rda" refers to 10 groups, 100 votes per session, 40 sessions
### correlations_hier.pdf and runtimes_hier.pdf are Figure 11


rm(list = ls())
setwd("~/Dropbox/hierIRT output")


runtimes <- matrix(0,nrow=3,ncol=4)
correl.all <- matrix(0,nrow=3,ncol=4)
colnames(runtimes) <- c("T=10","T=40","T=70","T=100")
colnames(correl.all) <- c("T=10","T=40","T=70","T=100")
rownames(runtimes) <- c("Court","Senate","House")
rownames(correl.all) <- c("Court","Senate","House")

load("G=10_J=100_T=10.rda")
runtimes[1,1] <- time[3]
correl.all[1,1] <- median(correl)

load("G=10_J=100_T=40.rda")
runtimes[1,2] <- time[3]
correl.all[1,2] <- median(correl)

load("G=10_J=100_T=70.rda")
runtimes[1,3] <- time[3]
correl.all[1,3] <- median(correl)

load("G=10_J=100_T=100.rda")
runtimes[1,4] <- time[3]
correl.all[1,4] <- median(correl)

load("G=100_J=500_T=10.rda")
runtimes[2,1] <- result[1]
correl.all[2,1] <- result[2]

load("G=100_J=500_T=40.rda")
runtimes[2,2] <- result[1]
correl.all[2,2] <- result[2]

load("G=100_J=500_T=70.rda")
runtimes[2,3] <- result[1]
correl.all[2,3] <- result[2]

load("G=100_J=500_T=100.rda")
runtimes[2,4] <- result[1]
correl.all[2,4] <- result[2]

load("G=500_J=1000_T=10.rda")
runtimes[3,1] <- result[1]
correl.all[3,1] <- result[2]

load("G=500_J=1000_T=40.rda")
runtimes[3,2] <- result[1]
correl.all[3,2] <- result[2]

load("G=500_J=1000_T=70.rda")
runtimes[3,3] <- result[1]
correl.all[3,3] <- result[2]

load("G=500_J=1000_T=100.rda")
runtimes[3,4] <- result[1]
correl.all[3,4] <- result[2]


## Bootstrap standard errors
load("G=10_J=100_T=10.rda")
sd(correl)
load("G=10_J=100_T=40.rda")
sd(correl)
load("G=10_J=100_T=70.rda")
sd(correl)
load("G=10_J=100_T=100.rda")
sd(correl)


runtimes <- runtimes/3600
runtimes
correl.all


pdf("correlations_hier.pdf")
t.vec <- c(10,40,70,100)
plot(0,0,type="n",xlim=c(10,100),ylim=c(0.8,1),
	xlab="Number of Sessions",
	ylab="Correlation", cex.lab=1.2)
points(t.vec, correl.all[1,], type="b", lwd=3, lty=1, pch=19, col="gray60")
points(t.vec, correl.all[2,], type="b", lwd=3, lty=2, pch=8)
points(t.vec, correl.all[3,], type="b", lwd=3, lty=3, pch=17)
legend(70,0.92, lwd=3, lty=c(3,2,1), pch=c(18,8,19), col=c("black","black","gray60"), bty="n", legend=c("G=500, J=1000","G=100, J=500","G=10, J=100"))
dev.off()



pdf("runtimes_hier.pdf")
plot(0,0,type="n",xlim=c(10,130),ylim=c(0,15), cex.lab=1.2,
	xlab="Number of Sessions",
	ylab="Computing Time (hours)")
points(t.vec, runtimes[1,], type="b", lwd=3, lty=1, pch=19, col="gray60")
points(t.vec, runtimes[2,], type="b", lwd=3, lty=2, pch=8)
points(t.vec, runtimes[3,], type="b", lwd=3, lty=3, pch=17)
text(120,13.8,"G=500, J=1000")
text(120,1,"G=100, J=500")
text(120,0.1,"G=10, J=100")
dev.off()


